Chimera patterns in conservative Hamiltonian systems and Bose–Einstein condensates of ultracold atoms

Experimental realizations of chimera patterns, characterized by coexisting regions of phase coherence and incoherence, have so far been achieved for non-conservative systems with dissipation and exclusively in classical settings. The possibility of observing chimera patterns in quantum systems has rarely been studied and it remains an open question if chimera patterns can exist in closed, or conservative quantum systems. Here, we tackle these challenges by first proposing a conservative Hamiltonian system with nonlocal hopping, where the energy is well-defined and conserved. We show explicitly that such a system can exhibit chimera patterns. Then we propose a physical mechanism for the nonlocal hopping by using an additional mediating channel. This leads us to propose a possible experimentally realizable quantum system based on a two-component Bose–Einstein condensate (BEC) with a spin-dependent optical lattice, where an untrapped component serves as the matter-wave mediating field. In this BEC system, nonlocal spatial hopping over tens of lattice sites can be achieved and simulations suggest that chimera patterns should be observable in certain parameter regimes.


Results
Nonlocal hopping model. Hamiltonian  where a i = √ n i e iθ i is a complex number representing the state of site i, such that |a i | is the amplitude, n i = |a i | 2 is the number of particles or density, and θ i is the phase. U is the nonlinear energy with the on-site nonlinear interaction U, and P is the hopping energy with the hopping strength P. G ij is the hopping kernel describing the G ij a * i a j Figure 1. Illustration of the dynamics of two different types of oscillators in a two-dimensional phase space. (a) A self-sustained oscillator with a limit cycle attractor. Trajectories near a limit cycle (represented by the dotted unit circle) move toward it. Energy disspation and driving are present such that two different initial states tend to the same asymptotic dynamics with the same oscillation frequency as time goes to infinity. (b) A conservative nonlinear oscillator. Typically, the oscillation frequency depends on the initial condition. Since energy is conserved, trajectories corresponding to different initial energies remain separated at all times. www.nature.com/scientificreports/ hopping from site r j to r i , with G ij = G ji . Typically, G ij decreases as the distance |r j − r i | increases and may be characterized by a hopping range R. For sufficiently small R, the hopping effectively becomes nearest neighbor. In this paper, we use G ij and R derived in Table 1. This Hamiltonian conserves both the energy and the particle number N = i n i . It can also be expressed using the canonical coordinate and momentum {q i , p i } , as well as action and angle variable {n i , θ i } (see SM Sect. S1). Note that the hopping term is quadratic a * i a j in the Hamiltonian, which is different from the usual quartic term of a particle-particle interaction n i n j for, say, the Coulomb interaction. Therefore, the corresponding dynamical equation contains the lowest order on-site nonlinearity and the nonlocal linear hopping term: where is the Planck constant, which we can set to = 1 without loss of generality by rescaling time. Note that this equation is the mean-field equation of the BHM with nonlocal hopping 55,56 . Moreover, the nearestneighbor variation of this equation is the discrete GPE 68 and the non-spatial variation is the discrete self-trapping equation 69 .
The dynamic equation of the NLHM can be rewritten in a dimensionless form using the rescaling a i → a i / √ n 0 , t → (Un 0 / )t , and P → P/(Un 0 ) where n 0 is the average number of particles per site. The equation becomes: iȧ i (t) = |a i | 2 a i − P j G ij a j , which depends only on the control parameters of rescaled hopping strength P and rescaled hopping radius R. Alternatively, Eq. (2) can be written in terms of θ i (t) and n i (t) as This explicitly shows that the evolution of the phase θ i (t) depends on the density n i (t) of the oscillators and vice versa. Even in the very weak hopping regime, they remain coupled to the lowest order. For dissipative systems as illustrated in Fig. 1a, if ṅ i ∼ 0 after dissipation in the weak coupling regime for all i, one can obtain a simplified phase dynamics. This is generally not possible for the conservative case with constant energy since, in general, a large n i at some site i has to be compensated by small n j at another site or sites to keep the energy constant. This highlights the important role of these conditions for conservative systems in contrast to dissipative systems. The dynamics of the NLHM can be found by solving Eq. (2) using standard numerical methods (see "Methods"), and the results for 1D and 2D are given in the following subsections.
Chimera patterns in 1D NLHM. An often used initial condition for chimera patterns is a random phase field 5,24,70,71 , in which chimera patterns can appear after a sufficiently long relaxation time. However, for NLHM, simulations show that the dynamics for such random initial conditions remains incoherent with no clear patterns over time. This is not unexpected since the spontaneous emergence of persistent patterns in spatially extended systems is typically tied to the notion of an attractor, which does not exist in our conservative model. Instead, incoherent and coherent regions-and, thus, chimera patterns-can sustain themselves over time as shown in Fig. 2a-d starting from initial conditions that are uniform with the exception of random phases (but not amplitudes or densities) in a small region. In particular, the time-averaged angular frequency �θ i � as shown in Fig. 2g is uniform in the coherent region and takes on a range of values in the incoherent region, thus, fulfilling the defining property of a chimera state. In terms of the temporal evolution, even though n i is constant initially, the random phases immediately induce fluctuations in the density as shown in Fig. 2b (see animations of the simulations in SM) as expected based on Eq. (3). Such a behavior can not be captured by simplified phase models by construction. To measure the coherence of the phase, we use the local order parameter O i = j G ij e iθ j

5
. The magnitude |O i | ∼ 1 when all phases θ j are the same within the hopping range R. As shown in Fig. 2e, |O i | takes on a minimum near the center of the incoherent region as expected. Moreover, the local order parameter does not converge but it keeps fluctuating as shown in Fig. 2f (see also Fig. 3f) due to the conservative nature of the system, which prevents relaxation behavior typical for dissipative systems. Figure 2a,c,e,g also indicate that the incoherent region is not fully desynchronized for this initial condition. A much stronger desynchronization can be obtained using an initial condition with both the phase and the amplitude random around the center region, see Fig. S2 in the SM. Hence, the initial amplitude plays a significant role for the characteristics of the observed chimera patterns in conservative systems, while this is typically not the case for dissipative systems, where initial fluctuations in the amplitude tend to be damped away. Another striking observation is that �θ i � can behave non-monotonically across the incoherent core (see Figs. 2g and also Fig. 3), whereas it typically changes Table 1. The D-dimensional hopping kernel G D (r) with r = |r j − r i | . K 0 is the modified Bessel function of the second kind. www.nature.com/scientificreports/ monotonically with distance from the incoherent core in the dissipative case and in simplified phase models 5 . The incoherent dynamics of the oscillators can be observed in Fig. 2h, where the trajectories of two oscillators inside the incoherent region are shown. The specific value of the hopping strenght P > 0 does not affect the chimera patterns qualitatively. However, for uniform initial conditions in the amplitude, the fluctuations in the amplitude can decrease when P decreases as shown in Fig. S3 in SM. While this could suggest that a simple phase description is sufficient in some special cases, such a simplification is generally not possible as already discussed above. Specifically, one distinctive feature of the NLHM is that the local phase oscillators can oscillate at any amplitude because of the lack of a limit cycle attractor. This can be observed using an initial condition with different amplitudes. An example is given in Fig. 3c,d where the initial density drops to zero and the phase changes by π at the center (this is a cross-section of a vortex phase  www.nature.com/scientificreports/ initial condition, see next section for more details). As suggested by previous studies 27, 28 , interesting chimera patterns can be formed spontaneously from such a regular initial condition. Here, the local phase incoherence and local density fluctuation around the center increase over time as shown in Fig. 3a,b. As Fig. 3h shows, the instantaneous frequency of the oscillators near the center can also change significantly over time. In particular, these oscillations have near zero amplitude as shown in Fig. 3b,h. In contrast, for the corresponding chimera patterns formed in dissipative systems with self-sustained oscillators the oscillations typically evolve close to the limit cycles in the weak coupling regime 6 .
Chimera patterns in 2D NLHM. Similarly to the 1D system, an initial condition with random phase regions can sustain itself over time in 2D. Here, we focus on such chimera patterns, in particular those where an incoherent region forms spontaneously around a phase singularity 25,27,28 . These patterns benefit from a topological protection in the sense that the incoherent core is robust against fluctuations in the phases. The first initial condition we examine is a spiral phase initial condition that is locally phase coherent everywhere except the center, with uniform density, as illustrated in Fig. 4a (see also "Methods"). With this initial condition, the system can spontaneously evolve into a state with a small incoherent core surrounded by a large spatially coherent region as shown www.nature.com/scientificreports/ in Fig. 4b for the phase field. Moreover, the density is randomized near the same core region in Fig. 4c. As shown by the dynamics of a cross-section in Fig. 4d, this spatial structure is sustained over long times (see Fig. S4 for snapshots and animations in SM). In addition, the same patterns can be observed even when the system size L and also R are increased (see Fig. S5 in SM). The local dynamics of the two oscillators in Figs. 4f,g clearly shows the difference between two regions: a i oscillates regularly far from the core, but not close to it. As in the 1D system, the incoherent region can only appear if the hopping range R is sufficiently large, here R 3 . Moreover, with nearest-neighbor hopping, the system reduces to the discrete GPE so that the incoherent region spreads out and interferes like a wave (see Fig. S8 in SM). All of these features are consistent with previous observations of chimera cores for driven-dissipative systems with self-sustained oscillators 24,25,28 . The distinct features in 2D are similar to the ones in 1D. This is shown in Fig. 4e-g for the angular frequencies and the trajectories in phase space. Note especially the strong variations in the average local rotation speed. In particular, the oscillators can exhibit significant variations in amplitudes as follows from Fig. 4c,g and the phase portrait in Fig. 4i that shows the phase and amplitude of all oscillators at a given moment in time. We would like to point out that after the formation of the chimera core, the pattern persists over the longest time scales we were able to simulate ( > 1000 spiral rotations). This observation suggests that if a random phase core is used as an initial condition, the chimera core pattern also persists over such long times scale. This is indeed what we observe (see Fig. S6 in SM).  www.nature.com/scientificreports/ The important amplitude-dependent dynamics without limit cycles can be clearly observed for the vortex phase initial condition with amplitude going to zero at center in Fig. 5 (see Fig. S7 for snapshots in SM), with a weaker hopping P = 0.1 . Similar to the 1D case discussed above, the fluctuations in the amplitude remain close to the initial condition for small P. In particular, oscillators with different amplitudes have different oscillating frequencies even in the weak hopping regime according to Eq. (3) with small corrections arising from the weak hopping. More importantly, as a conservative Hamiltonian system, it has time reversal symmetry and it conserves both quantities H and N (see Fig. S9 and animations in SM). This leads to persistent fluctuations or ripples as observed in Fig. 4b-d, which would be damped away in a dissipative system quickly. In addition, the results of the backward time evolution of the core region are very delicate. With a small perturbation, the background can evolve back to nearly the same state at t = 0 , but the core remains incoherent (see Fig. S9 in SM), which again signifies the difference between the two regions (see Sect. S4 and animations in SM). This suggests that the Poincaré recurrence time to a regular spiral-the time it takes to return within an arbitrarily small but finite distance to the original state (modulo possible rotations or translations)-is large and that the probability to encounter a regular spiral is zero in the infinite system size limit.
Moreover, the hopping energy P is not constant as shown in Fig. 4h even though the total energy H is constant. Hence, there is a conversion between P and U over time. This is different from a simple coherent and uniform distribution a i = √ n 0 having an energy per particle given by H/N = Un 2 0 /2 − Pn 0 with constant P and U . Note that all chimera patterns considered here do not correspond to ground states of the Hamiltonian but are excited states.
In realistic experimental systems, a small amount of particle loss typically exists and can be modeled phenomenonlogically by the term U → U − iU loss . Intuitively, the dynamics should not change significantly if the loss of the particles is less than half of the initial number of particles given by the condition U loss t/ 1 , Indeed, chimera www.nature.com/scientificreports/ patterns can, for example, still be observed with U loss /U = 0.02 at a sufficiently short time (Fig. S11 in SM).
Further details about such loss in 1D (Fig. S10 in SM) and 2D (Fig. S11 in SM) are discussed in Sect. S4 in SM.
Mechanism for nonlocal hopping and the minimal model. Mediating mechanism. The key idea for the mediating mechanism is to attach an inter-convertible mediating channel (labelled by ψ 2 ) to trapped states (labelled by ψ 1 ) as illustrated in Fig. 6a. With direct hopping, increasing the energy barrier between neighboring sites decreases both the hopping strength and the hopping range together. In contrast, if the particles can be converted into fast mediating states that do not experience any energy barrier, then the particles can physically jump much further away. Mathematically, this channel can be eliminated adiabatically (as done, for example, in 5,72 for non-Hamiltonian systems), resulting in an effective nonlocal model (see Fig. 6b) with independently adjustable on-site nonlinearity, hopping strength, and hopping range that can be tuned from nearestneighbor to global hopping.
Minimal model. A minimal mathematical model that captures the concepts of the mediating channel discussed above takes the form: for the localized ψ 1 and mediating ψ 2 components respectively. The corresponding Hamiltonian is given in Eq. (6) with appropriate parameters. The inter-conversion is governed by a detuning , which conserves the number of particles, and a coherent coupling with coherent oscillation frequency . This coupling may alternatively be referred to as Rabi coupling or Josephson coupling, depending on the physical systems being studied 62,63,73 . Eq. (4b) is essentially the Schrödinger equation for free particles with inverse mass κ = /(2m) > 0 and so the particles can propagate outward. The additional detuning in the far-detuned regime | | ≫ | | can ensure the mediating idea is well-defined: The number of particles N j = dr|ψ j | 2 in the mediating channel N 2 ≪ N 1 ≈ N can be neglected. Note that this model is not captured by the framework of nonlocal diffusive coupling 58 . It is (4a) i ψ 1 (r, t) =U|ψ 1 | 2 ψ 1 + �ψ 2 (4b) i ψ 2 (r, t) = − κ∇ 2 ψ 2 + �ψ 1 + �ψ 2 www.nature.com/scientificreports/ explicitly constructed to always preserve the conservation properties of the underlying Hamiltonian system, even when adiabatic elimination is applied.

Mediated hopping in ultracold atomic systems. Hamiltonian and dynamic equation. An ultracold atomic system of a general two-component GPE in a spin-dependent trap with coherent conversion is given by the Hamiltonian:
with and with the normalization N = N 1 + N 2 where N i = dr|ψ i (r)| 2 is the number of particles for each component. m i is the mass of the particles, V i (r) is the trap potential, g ij is the two-particle collision coefficient, and we assume g 12 = g 22 = 0 for the moment (see explanation below for non-zero case). The coherent oscillation term R represents the inter-conversion between the two components with the spatially homogeneous coherent oscillation frequency and the detuning i . By setting V i = 0 , m 1 → ∞ , and 1 = 0 , we arrive at the Hamiltonian for the minimal model discussed above. When a small nonlinearity exists in the mediating channel, the effective detuning becomes � → � + g 12 |ψ 1 | 2 + g 22 |ψ 2 | 2 if ψ i is uniform. Hence, the hopping radius decreases for g ij > 0 which is typical for atomic systems. Note that when |ψ i | 2 is small, the nonlinear effect can be ignored. It can be achieved by decreasing the density, which is one of the main technique used in the analysis of real systems below.
Mathematically, Eq. (4) can be obtained by setting appropriate parameters for the system described by Eq. (6). In particular, the absence of kinetic energy term in Eq. (4a) requires m 1 → ∞ . However, the mass m of interconvertible atomic systems are the same, so m i = m . To circumvent this, we can increase the effective mass; for example, by placing the atoms in a periodic lattice. This can be achieved by additionally setting V 2 = 0 , V 1 to be periodic, and 1 = 0 . Then the dynamic equation becomes 75 : (10a) i ψ 1 (r, t) = − κ∇ 2 + V 1 + g 11 |ψ 1 | 2 ψ 1 + �ψ 2 Figure 7. Chimera patterns in the minimal model with the direct simulation using Eq. (4) at t = 100 similar to Fig. 4b,c. The setting is the same as in Fig. 4  www.nature.com/scientificreports/ Only the positive detuning � = � 2 − ǫ 1 / > 0 is considered here as illustrated in Fig. 6c.
Mapping to effective NLHM. Note that direct adiabatic elimination does not work if states with high energy ǫ i>1 are occupied. This is because high energy states do not evolve slowly compared to the mediating component. To avoid occupying higher energy levels, we can confine the system to local ground states φ(r) with energy ǫ 1 and prevent excitation by choosing a suitable detuning such that ǫ 2 − ǫ 1 ≫ � ≫ |�| (see Fig. 6c). Under these constraints, along with adiabatic elimination, we can show (Sect. S2 in SM) that Eqs. (10a) and (10b) reduce to the exact form of Eq. (2) with U = g 11 |φ| 4 , P = � 2 /� , hopping kernel G D (r) in Table 1, and for d ≫ 2ℓ , where C D is a constant. Intuitively, particles staying in the mediating channel for a longer time have a larger hopping range R ∼ � −1/2 . Since the effective conversion region has a characteristic length scale 2ℓ in a unit lattice with length d, scaling with 2ℓ/d is expected. Indeed, we have the effective scaling � → � eff = (2ℓ/d) D � . The self-consistency condition for adiabatic elimination is ≫ Un 0 , P assuming all n i ∼ n 0 ( n 0 is the average number of particles per site). In this effective NLHM, a i in Eq. (1) represents the state of a localized wavepacket at site i. Moreover, the kernel G ij in Eq. (1) describes the matter-wave mediated hopping with wavepackets annihilated at site j and created at site i.
Optical lattice. The system discussed above requires a particle that is inter-convertible, which can be an atom with two different hyperfine states. A candidate is the Rubidium atom with hyperfine states |F = 1, m F = −1� and |F = 1, m F = 0� which has been realized in a spin-dependent trap 61 . Suppose the trapping potential is sinusoidal V 1 (r) = V 0 σ sin 2 (kx σ ) with wavelength , wavenumber k = 2π/ , lattice spacing d = /2 , and trap depth V 0 . The summation is taken over the lattice trap dimension as shown in Fig. 6c or e. For sufficiently large V 0 , all direct hopping can be suppressed, and the local ground states at trap minima can be approximated by a Gaussian φ σ (x σ ) = e −πx 2 /(2ℓ 2 σ ) / √ ℓ σ with ℓ σ = √ π /(mω σ ) . In this setting, the nonlinearity is enhanced by the high density since U = g 11 /W with effective volume W = 2 3/2 ℓ x ℓ y ℓ z . The constant can be found by numerical fitting, which gives C D ≈ 1 (see Sect. S3 and Fig. S1 in SM).
Achievable hopping range. For the hopping to be considered nonlocal, R > d must be satisfied. An example of Rubidium atoms is shown in Fig. 6d with d = 395 nm and a deep trap s = 40 (expressing V 0 = sE R in recoil energy E R = κk 2 ). With such a large s, as studied before 52 , the overlap between wavefuncion of neighboring cell is very small, the direct hopping is weak, and the system becomes a Mott insulator in the quantum regime. Nevertheless, mediated hopping can completely replace the direct hopping (with order R ∼ d , see Fig. 6d) and allow real time control. Since , , and U can be easily adjusted in experiments, there seems to be no upper bound on R. From a practical point of view, however, it is limited by the lifetime τ and experimental duration. A simple estimation of τ ∼ 1 s gives a maximum R ∼ 30d as shown in Fig. 6d.
Tuning nonlinearity and loss. The regime with competitive P ∼ Un 0 is the most interesting. However, a BEC in a 3D optical lattice using the parameters given above has a strong nonlinearity U/ = 2π × 2.23 kHz, which demands a large and, consequently, a small R. U can be reduced by the use of two tuning techniques: Decreasing the density, or utilizing the Feshbach resonance. The latter method can experimentally tune the nonlinearity over many orders of magnitude 66 . The former method is preferable because both nonlinearity and collision loss can be decreased simultaneously. In 1D and 2D lattices, the non-lattice dimension can be weakly trapped to reduce the density, resulting in a lattice of disk and cigarette-shaped wavefunctions respectively 76,77 . In this case, the dominant loss is the two-particle loss in the localized component. The rate of the two-particle loss can be estimated by U loss = L 11 /W and therefore half-life τ = W/L 11 with two-particle loss rate L 11 67 . This implies that τ ∼ ℓ z in 2D, so increasing ℓ z can improve the BEC lifetime.
Chimera patterns in BECs. The derivation of effective models implies that chimera patterns can also be observed in certain parameter regimes for Eqs. (4) and (10). The question is: can such parameter regimes be achieved in BEC experiments with current technology? The possible existence of chimera patterns in ultracold atoms is established in a parameter regime given in Fig. 8, based on a full simulation of Eqs. (10). Similar to Fig. 4, a random core appears eventually. Figure 8d shows the coherent oscillation between the two components with frequency ∼ √ 2 + 2 . Note that most of the atoms can be converted back after a full period, which confirms the physical picture discussed in Fig. 6a and is consistent with previous works 55,56 . The regime | | | | studied here is not in the far-detuned regime and may not be well described by NLHM, yet chimera patterns can still be observed in simulations. This suggested that chimera patterns do exist in a wide range of parameter regimes. As experimental techniques continue to improve, it will be possible to explore the adiabatic regime more closely.
Experimentally, the initial state can be prepared starting from a uniform BEC. Thousands of optical lattice sites [77][78][79] can be created with V 1 adiabatically turned on until the direct hopping is suppressed and mediated hopping begins to dominate. The energy shift induced by a short light-pulse can then be used to create any desired initial phase. The system states and dynamics may be detected by using various techniques such as optical readout, time of flight techniques, or matter-wave interference 80,81 . The loss U loss /U ≈ 0.017 here is comparable (10b) i ψ 2 (r, t) = − κ∇ 2 + � 2 ψ 2 + �ψ 1 www.nature.com/scientificreports/ with the discussion in the minimal model. Note that a small amount of loss can cause the BEC system to follow the classical trajectory 82 , and so each site can be well-described by a classical mean-field amplitude and phase. At the same time, our simulations suggest that chimera core patterns in 2D are particularly robust due to their topological structure. Specifically, if we start with a chimera core initial condition it can persist over long times. This is particularly useful if the lifetime of BECs is further limited in a given experiment by other experimental imperfections. All of this suggests that chimera patterns should be observable in experimental BECs.

Discussion
In summary, our work presents a Hamiltonian formulation of chimera states and demonstrates the existence of chimera patterns in three conservative Hamiltonian systems. The NLHM used in our study is a direct analogue of the nonlocal CGLE 5 , and our appoarch allows for the application of existing techniques in ultracold atoms and can be readily generalized to the quantum regime. Our simulations in realistic parameter regimes of BECs suggest that chimera patterns should be observable in experiments with ultracold systems using current technology. Additionally, our results suggest that the persistence of the incoherent region and the formation of a chimera core starting from a vortex or spiral initial condition in 2D are two distinct indicators of correct implementation of mediated nonlocal hopping, as opposed to local hopping which would result in the smoothing out of the incoherent region over time.
Our results in this paper are based on classical conservative Hamiltonian systems, which provides a new avenue to understand chimera patterns. These results may be extended into the quantum regime, since all of the physical processes that we analyzed are coherent and conserve both energy and particles. Equations (1)-(2), (4)-(10) can be quantized, and Eq. (1) becomes the Bose-Hubbard model with tunable mediated hopping 55 . This opens the door for the exploration of exotic condensed-matter states, such as supersolid states and quantum vortices with topological defect, in addition to other long-range effects 54,83 . The technique that we presented suggests that experimental studies of the synchronization and chimera patterns of a large number of oscillators may be feasible in quantum systems 37,[84][85][86][87][88] . We hope that our work here will motivate further studies of chimera patterns on ultracold atoms and quantum systems.

Methods
The numerical methods we used for solving Gross-Pitaevskii equations in Fig. 8 are the fourth-order time splitting method 89 . This method for the conservative systems automatically conserve the particle number. For all the other results, we used the standard fourth-order Runge-Kutta. The geometry used in the simulations is a square lattice with size L and the no-flux boundary condition. For the spiral initial condition, uniform density |a i | = √ n 0 is used and the state is given by a i (t = 0) = √ n 0 e i(k s r−tan −1 (y/x)) with r = (x 2 + y 2 ) . For the vortex-like initial condition, the state is given by a i (t = 0) = A i e i tan −1 (y/x) with A i = 1 − e −r/R vortex and R vortex is the length scale of the vortex. We use R vortex = R in this manuscript. For the system with a mediating channel, the channel is initially empty ψ 2 = 0.